program KineticQ
implicit none

character(len=79)ARGUMENT,RMSD
integer i,j,k,ios,nargs,totframe,NumA,NumO,NumH,prestate,PassTime2,PassTime3
integer AOroute,TotalAO,AHroute,TotalAH,OHroute,TotalOH,OAroute,TotalOA,HAroute,TotalHA,HOroute,TotalHO,AOHroute,TotalAOH,RouteSum,OAHroute,TotalOAH
integer,parameter :: first=999999
real,parameter :: Rcut=1.0
real :: tmp(7),RMSDA(first),RMSDO(first),RMSDH(first),ATP(first),Qi(first),Qf(first)
integer :: state(first),basin3(3),basin2(2),frame2(2),frame3(3),TauAO(first),TauOA(first),TauAH(first),TauHA(first),TauOH(first),TauHO(first),TauAOH(first),TauOAH(first)

nargs=iargc()
   if(nargs==1)then
     call getarg(1,RMSD)
   else
     write(*,*)
     write(*,*) "Wrong Input!"
     write(*,*) "./KineticQ.x Qs65.dat"
     write(*,*) 
     goto 100
   endif

ios=0
i=0
open(1,file=RMSD,status="old")
do while(.true.)
  read(1,"(A)",iostat=ios)ARGUMENT
	if(ios/=0)exit
	i=i+1
	read(ARGUMENT,*) (tmp(k),k=1,6)
	Qi(i)=tmp(2)
	Qf(i)=tmp(3)
!	RMSDH(i)=tmp(4)
!	ATP(i)=tmp(5)
enddo
close(1)

 totframe=i
 write(*,*) "There are ",totframe," frames!"

NumA=0;NumO=0;NumH=0
prestate=0
basin3(1)=0;basin3(2)=0;basin3(3)=0;basin2(1)=0;basin2(2)=0
frame2(1)=0;frame2(2)=0;frame3(1)=0;frame3(2)=0;frame3(3)=0
PassTime2=0;PassTime3=0
AOroute=0;OAroute=0;AHroute=0;HAroute=0;OHroute=0;HOroute=0
TotalAO=0;TotalAH=0;TotalOH=0;TotalOA=0;TotalHA=0;TotalHO=0
!AOHroute=0;OAHroute=0
!TotalAOH=0;TotalOAH=0

open(2,file="OUTPUT",status="replace")
open(100,file="OUTPUT100",status="replace")
!open(3,file="misfolded.dat",status="replace")
!open(4,file="ui.dat",status="replace")
!open(5,file="iu.dat",status="replace")
!open(6,file="ub.dat",status="replace")
!open(7,file="ib.dat",status="replace")
!open(8,file="HA.dat",status="replace")
!open(9,file="HO.dat",status="replace")
!open(10,file="AOH.dat",status="replace")
!open(11,file="OAH.dat",status="replace")
loop1: do i=1,totframe
!	if(RMSDA(i)<=Rcut .or. RMSDO(i)<=Rcut .or. RMSDH(i)<=Rcut)then
		if(Qi(i)<0.5 .and. Qf(i)<0.6)then
		 state(i)=1!!!Astate
		 NumA=NumA+1
 goto 66
		endif
		if( Qi(i)>=0.5 .and. Qf(i)<0.65 )then
		 state(i)=2!!!Ostate
		 NumO=NumO+1
 goto 66
		endif
		if(Qi(i)>0.75 .and. Qf(i)>=0.65)then
		 state(i)=3!!!Hstate
		 NumH=NumH+1
		endif

66		if(state(i)==0)then
!		 write(*,*) i-1," frame is misfolded!"
!		 write(3,*) i-1," frame is misfolded!"
		endif
		if(state(i)/=prestate)then
!		basin3(1)=basin3(2)
!		basin3(2)=basin3(3)
!		basin3(3)=state(i);
		basin2(1)=basin2(2)
		basin2(2)=state(i);
		frame2(1)=frame2(2)
		frame2(2)=i
		PassTime2=frame2(2)-frame2(1);
!		frame3(1)=frame3(2)
!		frame3(2)=frame3(3)
!		frame3(3)=i
!		PassTime3=frame3(3)-frame3(1);
			if(basin2(1)==1 .and. basin2(2)==2)then
			 AOroute=AOroute+1
			 TauAO(AOroute)=PassTime2
			 TotalAO=TotalAO+TauAO(AOroute)
!			 write(4,*) AOroute,frame2(1),frame2(2),TauAO(AOroute)
			endif
			if(basin2(1)==2 .and. basin2(2)==1)then
			 OAroute=OAroute+1
			 TauOA(OAroute)=PassTime2
			 TotalOA=TotalOA+TauOA(OAroute)
!			 write(5,*) OAroute,frame2(1),frame2(2),TauOA(OAroute)
			endif
			if(basin2(1)==1 .and. basin2(2)==3)then
			 AHroute=AHroute+1
			 TauAH(AHroute)=PassTime2
			 TotalAH=TotalAH+TauAH(AHroute)
!			 write(6,*) AHroute,frame2(1),frame2(2),TauAH(AHroute)
			 write(100,"(I6)") i
			exit loop1
			endif
			if(basin2(1)==2 .and. basin2(2)==3)then
			 OHroute=OHroute+1
			 TauOH(OHroute)=PassTime2
			 TotalOH=TotalOH+TauOH(OHroute)
!			 write(7,*) OHroute,frame2(1),frame2(2),TauOH(OHroute)
			 write(100,"(I6)") i
			exit loop1
			endif
!			if(basin2(1)==2 .and. basin2(2)==1)then
!			 OAroute=OAroute+1
!			 TauOA(OAroute)=PassTime2
!			 TotalOA=TotalOA+TauOA(OAroute)
!!			 write(7,*) OAroute,frame2(1),frame2(2),TauOA(OAroute)
!			endif
!			if(basin2(1)==3 .and. basin2(2)==1)then
!			 HAroute=HAroute+1
!			 TauHA(HAroute)=PassTime2
!			 TotalHA=TotalHA+TauHA(HAroute)
!!			 write(8,*) HAroute,frame2(1),frame2(2),TauHA(HAroute)
!			endif
!			if(basin2(1)==3 .and. basin2(2)==2)then
!			 HOroute=HOroute+1
!			 TauHO(HOroute)=PassTime2
!			 TotalHO=TotalHO+TauHO(HOroute)
!			 write(9,*) HOroute,frame2(1),frame2(2),TauHO(HOroute)
!			endif
!			if(basin3(1)==1 .and. basin3(2)==2 .and. basin3(3)==3)then
!			 AOHroute=AOHroute+1
!			 TauAOH(AOHroute)=PassTime3
!			 TotalAOH=TotalAOH+TauAOH(AOHroute)
!			 write(10,*) AOHroute,frame3(1),frame3(3),TauAOH(AOHroute)
!			endif
!			if(basin3(1)==2 .and. basin3(2)==1 .and. basin3(3)==3)then
!			 OAHroute=OAHroute+1
!			 TauOAH(OAHroute)=PassTime3
!			 TotalOAH=TotalOAH+TauOAH(OAHroute)
!			 write(11,*) OAHroute,frame3(1),frame3(3),TauOAH(OAHroute)
!			endif
!write(4,*) PassTime2,PassTime3,basin2(1),basin2(2)
		prestate=state(i)
		endif
!	endif
!  write(100,*) state(i)
enddo loop1
! write(2,*) "Basin Number : ","Apo: ",NumA,"Open: ",NumO,"Holo: ",NumH
! write(2,*) "Basin Probability :","Apo:",1.0*NumA/totframe,"Open:",1.0*NumO/totframe,"Holo:",1.0*NumH/totframe
! write(2,*) " "
! write(2,*) "Total Time: AOroute  OAroute  AHroute    HAroute    OHroute    HOroute"
! write(2,*) TotalAO,TotalOA,TotalAH,TotalHA,TotalOH,TotalHO
! write(2,*) " "
! write(2,*) "Number of AOroute  OAroute  AHroute    HAroute    OHroute    HOroute"
! write(2,*) AOroute,OAroute,AHroute,HAroute,OHroute,HOroute
! write(2,*) " "
! write(2,*) "Rate Constant KAB=NAB/TA WangYong-JCTC-ADK."
! write(2,*) "KAO  KOA  KAH  KHA  KOH  KHO"
! write(2,*) 1.0*AOroute/NumA,1.0*OAroute/NumO,1.0*AHroute/NumA,1.0*HAroute/NumH,1.0*OHroute/NumO,1.0*HOroute/NumH
! write(2,*) " "
! write(2,*) "Mean First Passage Time  MFPTAB=TotalAB/NAB"
! write(2,*) "AOroute  OAroute  AHroute    HAroute    OHroute   HOroute    AOHroute  OAHroute"
! write(2,*) 1.0*TotalAO/AOroute,1.0*TotalOA/OAroute,1.0*TotalAH/AHroute,1.0*TotalHA/HAroute,1.0*TotalOH/OHroute,1.0*TotalHO/HOroute,1.0*TotalAOH/AOHroute,1.0*TotalOAH/OAHroute
! write(2,*) " "
!	write(2,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
!	write(2,*) "ROUTE    COUNTS    AVERAGETIME(ps)"
22	if(AOroute>=0)then
	write(2,*) "u -> i: ",AOroute,1.0*TotalAO/AOroute
!	write(2,*) AOroute,1.0*TotalAO/AOroute
	endif
	if(OAroute>=0)then
	write(2,*) "i -> u: ",OAroute,1.0*TotalOA/OAroute
!	write(2,*) OAroute,1.0*TotalOA/OAroute
	endif
	if(AHroute>=0)then
	write(2,*) "u -> b: ",AHroute,1.0*TotalAH/AHroute
	endif
!	if(HAroute>0)then
!	write(2,*) "b -> u: ",HAroute,1.0*TotalHA/HAroute
!	endif
	if(OHroute>=0)then
	write(2,*) "i -> b: ",OHroute,1.0*TotalOH/OHroute
	endif
!	if(HOroute>0)then
!	write(2,*) "b -> i: ",HOroute,1.0*TotalHO/HOroute
!	endif
!!	if(AOHroute>0)then
!	write(2,*) "A -> O -> H: ",AOHroute,1.0*TotalAOH/AOHroute
!	endif
!	if(OAHroute>0)then
!	write(2,*) "O -> A -> H: ",OAHroute,1.0*TotalOAH/OAHroute
!	endif
	write(2,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"

!	RouteSum=AOroute+OAroute+AHroute+HAroute+OHroute+HOroute
!	write(2,*) "Fractional Flux  FAB=NAB/RouteSum"
!	write(2,*) "FAO  FOA  FAH  FHA  FOH  FHO"
!	write(2,*) 1.0*AOroute/RouteSum,1.0*OAroute/RouteSum,1.0*AHroute/RouteSum,1.0*HAroute/RouteSum,1.0*OHroute/RouteSum,1.0*HOroute/RouteSum
!	write(2,*) " "

close(2)
!close(3)
!close(4)
!close(5)
!close(6)
!close(7)
!close(8)
!close(9)
!close(10)
close(100)

100 end program
